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Abstract. The fast multipole method (FMM) performs fast approximate kernel summation to a specified 
tolerance e by using a hierarchical division of the domain, which groups source and receiver points into regions that 
' satisfy local separation and the well-separated pair decomposition properties. While square tilings and quadtrees 

^-H , are commonly used in 2D, we investigate alternative tilings and associated spatial data structures: regular hexagons 

■ (septree) and triangles (triangle-quadtree). We show that both structures satisfy separation properties for the FMM 
CN ■ and prove their theoretical error bounds and computational costs. Empirical runtime and error analysis of our 

^ , implementations are provided. 
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^ . 1. Introduction. In the original FMM works by Greengard and Rokhlin [7], multipole ex- 

^ \ pansions are executed along hierarchical centers of hypercube arrangements that span an input 

O ' domain. In the two dimensional case, the center of expansions can be mapped to a set of lat- 

tice points where their Vornonoi diagrams [Ij represent cell boundaries. That is, cell boundaries 
represent the equidistant points between nearest lattice points. For a square lattice configuration 
on the Euclidean plane, this equates to a square tiling and a general hypercube arrangement in 
^ ■ higher dimensions. Composing or decomposing these tilings naturally leads to a self-replicating 

hierarchical quadtree data structure [TT] that preserves local separation and well-separateness pair 
decomposition (WSPD) properties in [4j. This is necessary for the FMM to load-balance across 

■ levels and achieve linear runtime. 

o ■ 

^ ■ While the square lattice has been used by the FMM and the closely related treecode algorithm 

\ [2] since their inceptions, we are not aware of other lattice configuration having been explored. 

t> I In this paper, we investigate two alternative lattice groups that may be adapted into self-similar 

' arrangements; Regular triangular and honeycomb lattice groups translate to regular hexagonal and 

, triangle tilings. Although hexagonal tilings have been traditionally used in image processing jlOj 

and triangular tilings for mesh navigation [9J, they have not been applied to the FMM domain. We 
show how these arrangements form the bases for geometry preserving septree and triangle quadtree 
data structures that we adopt and implement for the hierarchical 2D FMM. Computational costs 
and error bounds induced by the data structures are derived to be comparable to that of quadtree. 

The outline of this paper is as follows: Section [2] provides an algorithmic preface of a hierarchical 
FMM for two dimensional coulombic equations. Section [3] illustrates separation properties that 
emerge from the tilings. Section [4] details a set of auxiliary functions shared amongst all data 
structures. Sections [5] and [6] introduce respective septree and triangle quadtree data structures and 
their implementations. Section [71 derives level invariant separation ratios. Section [8] provides an 
analysis of error bounds and computational costs. Section [9] presents empirical results that validate 
the theoretical analysis. Section [TOl concludes the paper and remarks on the generalizability of the 
data structures to higher dimensions. 
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2. Fast Multipole Method. The basic goal of the FMM is to evaluate the effects of a set 
of potentials denoted by source points Xi with strengths Ui on a set of receiver or target points yj. 
A potential is approximated upto an error bound via a series of hierarchical multipole expansions 
arranged over the spatial domain. In section [9l the algorithm is demonstrated on the 2D coulombic 
potential with kernel function 

^ji = log{yj-Xi), (2.1) 

over the complex plane though the discussion applies to general kernels. The total evaluation at a 
target point j is 

N 

$j = j = {1, . . . , m}, y.xeC. (2.2) 

i=l 

An overview of the FMM algorithm is provided in [3J. For brevity, the expansion and translation 
operators in [12j and [6J for the coulombic kernel are omitted. 

3. Separation Properties. A local separation property for multipole expansions and trans- 
lations guarantees a minimum separation distance between points assigned to non-adjacent tiles. 
Geometrically, two circles that circumscribe and inscribe a tile and its adjacent neighbors in Figs. 
13.11 bound a region. Formally, the separation ratio r/R oi radii between minor and major circles 
induces an error in the series expansion and so affects the number of truncation terms for a lower 
error bound provided in section [8l 




(a) Square (b) Hexagon (c) Triangle 



Fig. 3.1: Minor and major radii r and major R for the local separation property affects error in 
multipole expansions and translations 

The WSPD property defines a minimum separation distance between multipole to local (M2L) 
centers for translation. Formally, the distance between separated sets in Fig. 13.21 is expressed in 
terms of the ratio p/r and used to estimate the number of truncation terms. 




Fig. 3.2: WSPD distance p between lattice points affects error in M2L translations 
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4. Data Structures. To extend the local separation and WSPD properties beyond a basic 
tile, a concept of a cell is defined as a superset of tiles that shares a number of geometric and 
functional properties. For square and regular triangle tiles, cells are decomposable into smaller 
self-similar units. For square and regular hexagon tiles, cells are composable into larger aggregates. 
These hierarchical organizations give rise to indexable quadtree, septree, and triangle-quadtree data 
structures that when adapted for the FMM share the following attributes: 

1. Separation: Cells satisfy local separation and WSPD properties across all levels. 

2. Indexing: Cell indices satisfy some order relation in memory. 

3. Spatial addressing: Cell centers are computable from addresses and query points can find 
its bounding cell. 

4. Hierarchical addressing: Cell children, parent, and vertex neighbor indices are computable 
from addresses. 

The separation properties for septree and triangle-quadtree are derived in section [71 The in- 
dexing property for the quadtree and triangle-quadtree follow Morton Z curves in Figs. I4.1al and 
I4.1b[ The septree follows a spiral pattern in Fig. IS.lal The spatial and hierarchical addressing are 




(a) (b) 
Fig. 4.1: Morton Z curve indexing for (a) quadtree and (b) triangle quadtree 



accounted for by the following auxiliary functions: 

CellCenter(n, /). Returns the point coordinates of the center of cell n on level /. 
Celllndex(x, /). Returns a cell n on level / that contains point (x^y). 
Children(n). Returns cell indices of children cells of cell n. 
Parent (n). Returns a cell index of the parent of cell n. 

Neighbors(n, /). Returns cell indices on level / that share a vertex with cell n on level /. 
NeighborsE4(n, /). Finds cell indices required for M2L translation via 

Children[Parent{n) U Neighbor s{Parent{n), / — 1)] H Neighbors{n, I). (4.1) 

5. Septree. The septree data structure is an upward composable hexagonal tessellation of the 
Euclidean plane. While the basic hexagon units cannot be subdivided into self-similar components, 
they may be aggregated into larger hexagonal-like groups. A base 7 indexing scheme begins on 
level Ijnax with hexagon cell index O7 centered at the origin and adjacent to neighbors with indices 
{1, . . . , 6} tiled in a counter clock- wise direction. Aggregates of seven hexagonal cells on level Ij^ax 
form a single cell on level Imax — 1 where they can be tiled and indexed in a similar fashion as seen 
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in Fig. IS.lal Radial symmetry properties of the hexagon tessehation on level Imax give rise to the 
base 7 Generalized Balanced Ternary (GBT) [5j for cell indexing. This is investigated in section 
15.11 and its properties are exploited for a number of the FMM functions. 




(a) (b) 



Fig. 5.1: (a) Septree aggregates of hexagonal-like groups form a composable hierarchy where sym- 
metries along the vi, vs axis promote a base 7 GBT addressing scheme, (b) Table of all-pairwise 
GBT unit summations are analogous to vector summations to respective cells about the origin. 



5.1. Generalized Balanced Ternary. The generalization of Knuth's balanced ternary nota- 
tion hierarchically describes permutohedral regions of A^— dimensional spaces. In 2D, GBT indexing 
and arithmetic apply to hexagonal-like cells where the data can is decomposed. One important 
property realized by GBT addition is its duality with vector addition. That is, index arithmetic 
under GBT have a spatial analog with vector arithmetic in a coordinate system illustrated in Fig. 

EH 

For base 7 GBT unit addition, table ISTbl from Fig. IS.lal represents all pair- wise summations of 
unit vectors taken in each of the six neighboring cell directions. Subsequent summations of higher 
order cell indices are interpreted as a mapping between coordinates along axes vi^V2^vs and their 
base 7 cell addresses. For notation, denote index as the base 7 expansion of the usual base 10 
index n and the GBT addition operator as ®, e.g. the sum of vectors to cell indices 147 and SSr 
along V2 and vs axes is cell index 147 ® 357 — 327. 

Neighbors(n, /). The level independent neighbors of cell n are the pairwise GBT summations 
between cell index and the unit directions in Fig. I5.2a[ To adjust for the level, cell indices 
greater than 7^ — 1 are removed. The Neighbors(n, /) procedure is as follows: 

1. Let 7V7 = n7®{l,...,6} 

2. Return indices in A^io that are less than 7^ 

CellCenter(n, /). The center of cell n on level / is geometrically found via a series of transla- 
tions along zero-extended components of index n^. For notation, let index t — and component 
ti be the left-most digit. The cell center is expressed as 

1^1 

(x,^) = 5^r(t,z), (5.1) 

1=1 
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Fig. 5.2: (a) Neighbors of cell index 527 on level 2 are {65, 64, 6, 53, 50, 51}7. (b) Cell center of 
index 327 found via translations along zero-extended components 3O2 and 27. Zero-extended index 
6O7 found via 67 ® I7 ® I7 



where \t\ is total number of components and function T : ^ extends the i^^ component of 
t with \t\ — i zeros and returns the corresponding vector to the cell from the origin, e.g. the cell 
center of index 3427 is the sum of vectors to the cell centers of indices {300, 40, 2)7. 

The zero-extended i^^ component is the GBT summation of three shorter zero-extended com- 
ponents. The base case in table ISTSal shows how successive GBT summations of a unit index j and 
twice 1 + (j mod 6) yield a zero-extended jO. For z < zero-extending the components ignores 
the preceding zeros, e.g. IO7 ® 2O7 ® 2O7 = IOO7. This is expressed by the relation 

t,0\t\-' = t.ol'l-^-^ ® [1 + {U mod 6)]0l'l-^-^ ® [1 + {U mod 6)]0l'l-^-^ (5.2) 



1®2®2 = 10 
2 ® 3 ® 3 = 20 
3®4® 4 ^30 

4 ® 5 ® 5 = 40 

5 ® 6 ® 6 = 50 
6® 1® 1 ^60 

(a) (b) 

Fig. 5.3: (a) GBT unit addition for zero-extension, (b) Geometric analog for a zero-extension via 
translation b 




A geometric interpretation of index t^O'^' ^ in Fig. 15. 3b I reveals that the translation b is the sum 
of two vectors separated by 27r/3 radians with a two-to-one magnitude ratio. Solving for the angle 
of separation and the magnitude of vector b in appendix eqs. ID. II and ID. 21 yield 9 = arctan ^ 
and b = a\/7. Hence, successive zero-extensions of a component ti are rotated by 9 radians with 
magnitude a^+i = a^\/7 with a base magnitude ao = 
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The function T{t, i) can now be expressed in polar coordinates 

j = {\t\ + [Imax - I), 



\/3 vr TT 

7 arctan (- ti , Ri 

J 2 3 6 



(5.3) 



before T{t^i) — {Ri cos 0i, Ri sin. 0i) gives the translation in a Cartesian coordinate system. The 
CenCenter(n, /) procedure is as fohows: 

1. Let = (0,0), i = 1 

2. Update (x, y) = (x, y) + {Ri cos Ri sin^^) from eq. 15.31 

3. Increment i and repeat from step [2] until i = / 

4. Return point (x, 

Celllndex(x, y^ I). An approximation of the cell that contains the query point (x, y) is found 
via a change of basis 

/ 2x y\/3 — X \ 
3^' 3r ' 



{b,c) 



(5.4) 



where coordinates (6, c) are defined along axes V2,V3- These associated cells along the axes are 
erroneous as points near the corners of true bounding hexagon may project incorrectly in Fig. I5.4al 




Fig. 5.4: (a) Transforming query point (x, y) (6, c) via change of basis may lead to incorrect cells 
along V2^vs axes, (b) Hexagonal edges are Voronoi diagrams of equilateral triangular lattice points 
that construct cell centers 



To address this issue, boundary conditions are checked by observing that hexagonal edges are 
the Voronoi diagrams of equilateral triangular lattice points. That is, hexagonal edges represent 
points equidistant between nearest lattice points. The query point (x, y) that falls within a hexagon 
would be the nearest neighbor of its center or lattice point. Furthermore, the query point may be 
bound between four neighbouring candidate lattice points in Fig. I5.4b[ If the lattice points are 
defined w.r.t. axes V2^vs as 

(Srb rV3{2c + b) \ 

{X,y)lattice = I ^ 1 , (5.5) 
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then the candidate lattice points have coordinates 

{mAc\),m,[c\),m,\c]),{\b-\,\c\)}. (5.6) 

Finding the nearest lattice point w.r.t. the query point (6, c) after converting back to Carte- 
sian coordinates via eq. 15.41 yields the hexagon and its respective cell index t — nj. The first / 
components or ti-i is the the cell index on level /. The Celllndex(x, /) procedure is as follows: 

1. Compute axes coordinates (6, c) from eq. 15.41 

2. Compute z = 1 : 4 candidates lattice axes coordinates (6, c)i from eq. 15.61 

3. Convert candidates lattice coordinates to Cartesian coordinates (x, y)i with eq. 15.41 

4. Let the nearest candidate lattice point w.r.t. query point {x^y) be (x,^) and equivalently 
(6,c) 

5. Find cell indices u^^ for coordinates 6, c along axes vs 

6. Return u® v 

6. Triangle Quadtree. A variation of the original quadtree data structure consists of a regu- 
lar triangular tessellation of the Euclidean plane that is strictly decomposable. While each triangle 
may be subdivided into four similar units (center, vertical, left, right), a group of triangles is not 
upward composable in the sense that vertical orientations of descendant triangles depend on the 
orientation of their common ancestor triangle or root. That is, all children of a center cell will 
have inverted its orientation. For notation, denote the default upright orientation as upi = 1 and 
the inverted orientation as upi = —1 for a triangle on level i. Triangle quadtree indexing begins 
on level i = with cell index centered at the origin and enclosing the entire domain. Children 
cells on level z + 1 with indices {0, . . . , 3} are associated with cell types {center, vertical, left, right} 
ordered w.r.t. the parent cell. 

Neighbors(n, /). For regular quadtrees, neighbors of a cell index n need only to share an 
edge. For triangle quadtrees, denote the left, right, and vertical adjacent neighbors as {L, R, V} in 
Fig. [6Tal 




(a) (b) 

Fig. 6.1: (a) Neighbor finding from cell index 03 in the vertical direction. Cell is nearest ancestor 
with vertical sibling cell 1 on level 1. Reflect the path leading to ancestor cell gives cell index 13. 
(b) Neighbor finding from recursive calls to adjacent neighbors 

Similar to [9j, the adjacent neighbor of cell n in direction k is found with an ancestor-sibling- 
reflect method in the adjacentNeighbor(n, k) procedure: 

1. Search for the closest common ancestor containing a sibling in direction k. For notation, 
let cell index t = 714 and matrix D in table 16.11 map cell source types and search directions 
to destination cell types. From right to left components find the first valid entry Dt.^j^. 
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Source cell type 


Left cell ID 


Right cell ID 


Vert cell ID 





2* 


3* 


1* 


1 


3 


2 


0* 


2 


1 


0* 


2 


3 


0* 


1 


3 



Table 6.1: Mapping from source cell types and search directions to destination cell types. * entries 
indicate valid siblings of a source cell 



2. Move to sibling node. Replace component ti with entry Df.^j^. 

3. Reflect the path taken to reach ancestor. For j = {(z + 1), . . . , Imax}^ replace component tj 
with entry Dt-^k- 

4. Return cell index t 

To find all neighbors that share a vertex with cell n, make a set of recursive calls to adjacentNeighbor(n, 
k) in Fig. I6.1b[ The Neighbors(n, /) procedure is as follows: 

1. Let {L, R, V} be adjacent neighbors of cell n in directions left, right, vertical 

2. Let {VL^VR, LL, RR, LV, RV} be adjacent neighbors of cells V,L,R in directions left, 
right, vertical 

3. Let {VLL^VRR, LVR} be adjacent neighbors of cells VL^VR^ LV in directions left and 
right 

4. Return ceh indices {L, R, V, VL, VR, LL, RR, LV, RV, VLL, VRR, LVR} 

CellCenter(n, /). The center of cell n on level / is found via a series of translations for each 
component of index 714. Let index t = 714, component ti be the left-most digit, and cell center 

1^1 

= (6.1) 

i=l 

where function F : Z'^ ^M? returns a vector to cell index U that is adjusted for level and orientation, 
e.g. the cell center of index 3024 is the sum of vectors to cell centers of indices {3, 0, 2)4 on levels 
{1, 2, 3} with orientations {1, —1, —1}. 

To determine the orientation of index t on level /, observe that the orientation only flips for 
center cells in Fig. I6.2a[ This is the zero-parity of cell index t (even is 0, odd is 1) adjusted for 
level / and the number of components \t\. The function 

isUp(t, /) = 2[(parity(t) + / - |t| + 1) mod 2] - 1, (6.2) 

returns 1 for the up orientation, and —1 for the inverted orientation. To determine the orientation 
of component make a query to isUp(ti:^,/ — \t\ + i) where cell index ti:i is the level adjusted 
ancestor. 

Note that triangle centers for components U = are equivalent to its ancestor's center as no 
translations are necessary. The function F(t, i) can now be expressed in polar coordinates 

j = / - |t| + z, 1 < < 3, 

^. = isUp(t,.,,)(| + ^^^), R. = 2^—^r, ^'-'^ 

where radius R is scaled by a factor of two per level in Fig. I6.2a[ The CellCenter(n, /) procedure 
is as follows: 
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Fig. 6.2: (a) Triangle orientation determined by cell index zero-parity. Translation magnitudes to 
adjacent neighbors scale by a factor of two per level, (b) Triangle edges are Voronoi diagrams of 
hexagonal lattice points at cell centers. A cell index is found via searching nearest lattice points 
within a current triangle. 



1. Let current center coordinates (x^y) = (0,0), z = 1 

2. Update center (x, y) = (x, y) + (Ri cos 0i, Ri sin 9i) from eq. 16.31 

3. Increment i and repeat from step [2] until i — I 

4. Return (x, y) 

Celllndex(x, y^ I). A cell index t — nj contains point (x,y) on level / if all ancestor cell 
indices ti-i also contain it. Since triangles descendants are contained within ancestor triangles, a 
search proceeds from level to / s.t. successive descendant triangles contain point (x,y). 

Observe that triangle edges form Voronoi diagrams of regular hexagonal lattice points. That 
is, the triangle edges represent points equidistant between nearest lattice points. Hence, a query 
point (x, y) that falls within a triangle is also the nearest neighbor of the triangle's center or lattice 
point in Fig. I6.2b[ 

To track the progress of triangle centers, denote point {u^ v) as the current center. Note that 
triangle centers for components = are equivalent to its ancestor's center. That is, if a center 
triangle is found to be the nearest neighbor, then the current center remains unchanged while the 
vertical orientation up flips. In polar coordinates, denote the three triangle lattice points about the 
origin as 



where similar to the CellCenter function, the radius R is scaled by a factor of two per level. The 
Celllndex(x, /) procedure is as follows: 

1. Set to origin the current center (u^v) = (0,0), let z = 1 

2. Assign component U the cell type with the nearest lattice to query point (x, y) on level i 

3. Update current center (u^v) = (u^v) + (Ri cos Ri sin 9i) from eq. 16.41 

4. If ti = 0, then flip orientation upi = —upi-i 

5. Increment i and repeat from step [2] until i = / 

6. Return t 




— 2, 




(6.4) 
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7. Separation Ratios. To show that septree local separation and WSPD ratios are level 
invariant, consider the base case in Fig. 13. Ibl where ratio r /R— 1/2. For successive levels, lattice 
points induce Voronoi diagrams that form rotated and scaled hexagonal boundaries that intersect 
points along minor and major radii r and R in Fig. I7.1a[ The rotation and scaling are proportional 
to eq. 16.31 adjusted for level i. The WSPD distance p = r in Fig. 17. Ibl and local separation ratio 
are preserved as the induced hexagonal boundaries are self-similar to the base case. 




(a) (b) 



Fig. 7.1: (a) Local separation ratio r /R drawn from level 1 lattice points and their induced hexag- 
onal Voronoi diagrams. Either radius r or i? is computable from magnitude h obtained from eq. 
16. 3[ (b) WSPD distance p — ?>r — a — r generalized over level induced hexagonal boundaries 

To show that the triangle-quadtree local separation ratio is level invariant, observe that ratio 
r/i? = 1/2 for the same level in Fig. I6.2a[ On successive levels, a major radius R is equivalent to 
the minor radius r from the preceding level and can be computed from eq. 16. 4[ The WSPD ratio 
p/r — \/7 — 2 in Fig. 17.21 is derived in appendix [El Similar to the quadtree, the level invariance 
stems from self-similar tiling arrangements with the base case in Fig. I3.1cl 




(a) (b) 
Fig. 7.2: (a) Triangle tilings for nearest M2L translation, (b) Geometric analysis for WSPD distance 

8. Cost Analysis. To estimate the costs of a hierarchical FMM, a uniform source and target 
point distribution over the data structure's geometric domain is assumed [8j. The total computa- 
tional costs w.r.t. the quadtree, septree and triangle quadtree are derived in appendices lAj iBl and 
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struct 


P4 


P2 


r/R 


p/r 


Sopt 


Opt Costs 


Q-tree 


27 


9 




2{V2-1) 


Cii6ivy& 

\ 27M J 


{M + N + 22M{MN)-^)p 


S-tree 


42 


7 


1 
2 


1 


\ 42M ) 


{M + N + 35.2{MNy^)p 


TQ-tree 


39 


13 


1 

2 


V7-2 


\ 39M / 


{M + N + 4Q.Q5{MNy^)p 



Table 8.1: FMM properties for quadtree, septree and triangle quadtree where p is the number of 
truncation terms for each expansion, P4 is number of M2L translations per cell, P2 number of cells 
in neighborhood, p is the multipole to local separation distance and Sopt the optimal number of 
points per cell 



ICl The overall costs depend on a density quantity Sopt of source and target points per cell. This 
density quantity is a function of the number of M2L translations P4 and the number of neighbor 
cells P2 which is optimized by equating the number of M2L translations to the number of direct 
evaluations per cell. Substituting density Sopt back into the FMM total costs yields the optimal 
FMM costs hsted in table EH 

The total FMM costs also depend on the number of truncation terms p which are chosen 
a priori to guarantee a minimum error bound. This quantity depends on the kernel expansion 
and translation formulations as specified in [12j. Maximum absolute error bounds [6J are written in 
terms of local separation ratio r/R and WSPD ratio p/r between M2L translations. For a multipole 
expansion and translation, the absolute error is bounded by 

' ^' - R-r \rJ 



< Lh]Ll ( 1 . (8.2) 



For a M2L translation, the absolute error is bounded by 

TT^/r 

Last, the local to local translations are exact. 

For a fixed number of truncation terms, the total optimized cost for septree is slightly greater 
than that of quadtree as the lower neighbor count does not fully compensate for the greater number 
of M2L translations per cell. The triangle quadtree is less cost efficient in both categories and so 
yields a larger leading coefficient term. For multipole expansions, the larger local separation ratios 
r/R in both septree and triangle quadtree suggest wider error bounds than that of quadtree. For 
M2L translations, the smaller WSPD ratio p/r for septree suggests a tighter error bound than that 
of both quadtree and triangle-septree. This last property may have considerations in setting the 
maximum l^ax of the data structure. 

9. Experiments. To validate the theoretical FMM costs and error bounds from section [HI 
a set of test cases suitable for uniform conditions for each data structure is generated. That is, 
the input domain consist of Sopt uniformly random source and target points fitted to each cell. A 
regular polygon point picking method in Fig. 19. II handles the various tile geometries. For a regular 
cell, the polygon is subdivided into disjoint isosceles triangles per edge where halves of triangles are 
rearranged to form rectangles. Uniformly random points (x, y) are picked from this rectangle and 
a third uniformly random variable z assigns the point to one of the isosceles triangles. To modify 
the total number of source and target points, the total number of occupied cells is reduced. 
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Fig. 9.1: Regular polygon point picking via parameterization into uniform variables (x^y^z) where 
a point (x, y) is uniformly sampled from a rectangle while uniformly random variable z chooses the 
triangle 



For experiments, both runtime and error estimates of the FMM are considered while fixing either 
the number of source and target points {N, M), or the number of truncation terms p. Runtime 
comparisons w.r.t. the direct evaluation method are shown when appropriate. The implementation 
is written and tested on Matlab 2010b and Intel i7-2630QM hardware. 

The runtime results for a variable number of source and target points in Fig. 19.21 match 
the theoretical estimates. All three hierarchical FMM data structures obtain linear asymptotic 
runtimes. The basic quadtree data structure outperforms the septree by a small margin for both 
max levels Imax- The triangle quadtree performs the slowest out of the three. The direct method 
obtains a quadratic runtime and where an estimated cross-over point between the number of source 
and target points is between 10^*^ to 10^'^^. The error results in Fig. [931 indicate a linear increase in 
error w.r.t. the number of source and target points N and M. The quadtree obtains the least error 
presumably due to the smallest local separation ratio. One explanation for the greater maximum 
error in the septree compared to the triangle quadtree is the larger number of corner points in the 
tiling where source-target point distances are minimized. 

The runtime results for a variable number of truncation terms p in Fig. 19.41 reveal a cross- 
over point between septree and triangle-quadtree. This may be due to the greater number of 
M2L translations in the septree and the inefficiencies in the implementation where the translation 
matrices are computed on the fly. The error results in Fig. 19.51 exhibit an exponential error loss for 
increasing truncation terms. 

10. Conclusions. In this paper, we have shown that regular hexagonal and triangle tilings 
of the Euclidean plane generate septree and triangle-quadtree data structures that satisfy local 
separation and WSPD properties for the FMM. We derived their implementations and respective 
geometric properties with regards to FMM error bounds and costs. The empirical results validated 
the theoretical claims and have shown to be comparable to the original quadtree. While the quadtree 
remained ideal for a uniformly distributed data set, both the septree and triangle-quadtree may 
have applications for non-uniform domains, especially when the input is clustered about respective 
lattice points. 

As a final remark, the septree's base hexagon unit is an order-3 permutohedron. Permuto- 
hedrons are (d — 1) dimensional polytopes embedded in d dimensional space that can tessellate 
the domain and is indexable via a similar GBT addressing scheme. Future work may investigate 
constructions of a permutohedron-tree and its separation properties in higher dimension. 
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Runtimes for Fixed Truncation Terms p=12 
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Appendix A. FMM Quadtree Costs. FMM evaluation costs before and after substituting 
optimal point density per cell for quadtree data structure in eqs. lA.ll and IA.2[ 

P4(d) = 3^(2^-1), P2(d) = 3^. 



cost(FMM) = (M + N)P + I^K^;^iP,{d) + 2) ^ J P 

+ M{P2{d)s + P). 



opt cost = (M + N)P + \^ + y^i^ I {MNin"^^ - 1) + 2)3'')-^P 
[M + N + 32.64(MAr)-^) P, d = 2. 



(A.l) 



(A.2) 



Appendix B. FMM Septree Costs. FMM evaluation costs before and after substituting 
optimal point density per cell for septree data structure in eqs. IB.ll and IB.21 



N = 7^% S = 7^% L = L^-Ls, K = 7^, P4 = 42, P2 = 7. 

/ 7 74 _ p 72 \ 
cost(FMM) = (M + A^)P+ f K-(P4 + 2) j P^ + M(P2S + P). (B.l) 

opt cost = (M + N)P + l^yi + (308MAr)-5p 
« (M + AT + 35.2(MAr)-5) p 
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Fig. 9.4: Runtime (seconds) for FMM data structures with variable number of truncation terms p 
and fixed number source and target points N — M — 5120 



Appendix C. FMM Triangle Quadtree Costs. FMM evaluation costs before and after 
substituting optimal point density per cell for triangle quadtree data structure in eqs. IC.ll and lC.2[ 

AT = 4^*, 5 = 4^% L = L*-L„ K = 4^, P4 = 39, P2 = 13. 



cost (FMM) = (M + N)P + ( K-{P^ + 2) - 

3 3 



44 _ p,72\ 

\ p2^^^p^^^py ^c.l) 



opt cost = (M + N)P + ( Y ^ + Y ^ ) (533M7V)-^P 



(C.2) 



(M + TV + 46.65(MiV)-^) P. 



Appendix D. GBT Translations. Septree GBT translations for zero-extending indices. The 
angular separation 6 in Fig. I5.3bl is 



sin6> _ sin 120 _ sin (60 - 9) 
2a h a 



2 



(D.l) 



The magnitude of translation vector b in Fig. I5.3bl is 



2 



2 



(D.2) 
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Fig. 9.5: Maximum absolute errors for FMM data structures with variable number of truncation 
terms p and fixed number source and target points N — M — 5120 



Appendix E. Triangle-quadtree WSPD. The min M2L translation distance p for triangle 
quadtree in Fig. 17.21 is 

sin (9 sin 150 30 - 6> ^ \/3 

= = -=- ^ — arctan , 

r a r\/3 9 

r V3 — 1 

- = a sin arctan = a — ^ = = a — —, (E.l) 

2 9 /i+r^f 



V 9 ; 

a r\/7, p r{V7 — 2). 



